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The equilibrium and fluctuation methods for determining the surface tension, a, and bending 
modulus, k, of a bilayer membrane with a fixed projected area are discussed. In the fluctuation 
method the elastic coefficients a and k are measured from the amplitude of thermal fluctuations 
of the planar membrane, while in the equilibrium method the free energy required to deform the 
membrane is considered. The latter approach is used to derive new expressions for o and k (as well 
as for the saddle-splay modulus) , which relate them to the pair-interactions between the amphiphiles 
forming the membrane. We use linear response theory to argue that the two routes lead to similar 
values for o and ft. This argument is confirmed by Monte Carlo simulations of a model membrane 
whose elastic coefficients are calculated using both methods. 
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I. INTRODUCTION 



The Bilayer membrane, a double sheet of surfactants separating two aqueous phases, is one of the structures formed 
by the self-assembly of amphiphilic molecules in water . The driving force in this process is the hydrophobic effect 
which favors exposing the hydrophilic part of the molecules to the water while shielding the "oily" part from aqueous 
contact The ongoing interest in such membranes is due to many reasons, among which are their predominant 

role in the organization of the biological cells , and their various applications in many industrial sectors [j| . Bilayer 
amphiphilic sheets have very special mechanical properties: While being strongly resistant to lateral mechanical 
stretching or compression, they are highly flexible and can exhibit large thermally excited undulations 0,0- This 
unique elastic behavior, namely the stability against external perturbations on the one hand, but the case in going 
from one shape to another on the other hand, is important for the activity of living cells §]. Consequently, there has 
been a great effort to understand the elasticity of bilayer systems 0, 13 EH • 

Bilayer membranes are quasi two-dimensional (2D) objects: their thickness is typically of the size of a few nanometers 
(roughly, twice the length of the constituent amphiphilic molecules), while their lateral extension can reach up to 
several micrometers. Since the membrane appears as a thin film on the mesoscopic scale, its physical properties are 
often studied using coarse-grained phenomenological models treating the membrane as a smooth continuous 2D sheet 
0,0,0; El- Membrane elasticity has been traditionally studied using the Helfrich effective surface Hamiltonian which 
relates the elastic energy to the local principle curvatures of the membrane c\ and c 2 , and which has the following 
form [U: 



n= ds 

J A 



i«o (J - 2c ) 2 



k K 



(1) 



where J = C1 + C2 and K = c\Ci are the total and Gaussian curvatures respectively. The integration in Eq.JQl is carried 
over the whole surface of the membrane. The Helfrich Hamiltonian is derived by assuming that local curvatures are 
small, and the free energy can be expanded to second order in J and to first order in K. It, therefore, involves four 
phenomenological parameters: the spontaneous curvature Co, and three elastic coefficients - the surface tension oo, the 
bending modulus kq, and the saddle-splay modulus Rq, whose values depend on the area density of the amphiphiles. If 
the number of these is fixed, then one should also consider the corrections to Hamiltonian Q due to the changes in the 
area of the fluctuating membrane. For weakly fluctuating membranes these corrections can be assumed to be small. 
The surface tension, which is usually associated with the free energy cost for adding molecules to the membrane (at 
a fixed density), is related in the case of membranes with fixed number of amphiphiles to the area-density dependent 
(Schulman) elastic energy [l3L liH HH . 
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The Helfrich Hamiltonian has been very successful in describing the shape and the phase diagram of complex 
interfaces Il6l ll7lll8j. It also yields a correct description of the thermal fluctuations around the equilibrium surface 
state [TflL l2fj EH] . and of the entropic forces between membranes j^- Because it is phenomenological, the Helfrich 
Hamiltonian provides no information about the values of the elastic coefficients. Many theories have been developed 
that attempt to relate the elastic coefficients introduced by the Helfrich Hamiltonian to microscopic entities and the 
interactions between them pii 123 , EEL EfiLE^ . In fact, these theories are usually concerned with the free energy of the 
surface, rather than the Hamiltonian. The free energy is assumed to have the same form as the Helfrich Hamiltonian 
and, hence, usually called the Helfrich free energy (see a more detailed discussion in section ITTll . The coefficients 
appearing in the expression for the free energy, which we denote by a, n, and R, are also referred to as the surface 
tension, the bending modulus, and the saddle-splay modulus, respectively. Despite the similarity in names, there is 
a significant difference between the Hamiltonian coefficients (with the subscript 0) and the free energy coefficients. 
The former are "material properties" which depend on the internal (potential) energy of the surface. The latter, on 
the other hand, are thermodynamic quantities and, as such, are also influenced by the entropy associated with the 
thermal fluctuations of the system. Their values, therefore, may also depend on the temperature and the size of the 
system. 

In addition to the above mentioned theories, there has been also an effort to analyze the elastic behavior in the 
context of the thermodynamics and statistical mechanics of curved interfaces [28l I2H IsH . Ell E2I Isflj . The last 
approach has the potential of providing exact "virial" expressions for a, k, and R in terms of the microscopic forces 
between the amphiphiles and the pair distribution function. One of the systems whose statistical mechanics has been 
studied extensively is that of a simple liquid-vapor interface. Although this seems to be a rather simple system, the 
determination of its elastic moduli is quite complicated and involves a set of technical and conceptual problems. Below 
we discuss some of them: 

• One problem is related to the finite thickness of the interface, namely to the fact that the local concentration 
is not a step function but changes gradually while going from one phase to the other. Consequently, there 
is some ambiguity about the location of the dividing plane that separates the two phases and to which the 
Helfrich Hamiltonian is applied. It turns out that the values of the rigidity constants k and R (the coefficients 
of the second order terms in the curvatures c\ and c%) depend on the choice of the dividing surface 0. The 
dependence of the rigidity constants on the reference surface had led people to question the validity of continuing 
the Helfrich free energy expansion beyond the linear term in curvature. This problem has been recently tackled 
by van Giessen and Blokhuis |35j| who used computer simulation to determine the rigidity constants of a curved 
liquid-vapor interface in a system of particles interacting via a truncated Lennard-Jones (LJ) potential. They 
have demonstrated that although one needs to state which convention for locating the dividing surface is used 
when providing the values of k and R, this fact does not render the Helfrich free energy useless, nor does it 
diminish the importance of these quantities in describing the elastic properties of the interface. 

• A second problem that makes the determination of the rigidity constants difficult is a technical one: In their 
paper van Giessen and Blokhuis used the virial expressions given in Ref. |36| to evaluate the values of k and R. 
These expressions relate the rigidity constants to the derivative of the pair density distribution function with 
respect to the radius of curvature R c . This means that the values of the rigidity constants of a planar interface 
cannot be determined from the simulation of that system only, but it is necessary to perform a set of simulations 
of curved interfaces with very large values of R c . For the interfaces investigated in Ref. [35|. it turns out that 
in the large R c regime the dependence of the pair density function on R c is very weak. Consequently, it was 
impossible to determine k and R accurately, and only a rough estimate of these quantities could be obtained. 

• A third problem, a more fundamental one, is related to the method of calculating the rigidity constants k 
and R and to our interpretation of their physical meaning. The theoretical and experimental methods for 
determining the elastic coefficients of interfaces can be classified into equilibrium (or mechanical) methods and 
fluctuation methods |37l l38|| . The difference between these two approaches is in the context in which the 
Helfrich Hamiltonian and the associated free energy are used: In the equilibrium approach one extracts the 
elastic coefficients by comparing the free energies of two equilibrium surfaces with different curvatures. In the 
fluctuation approach, on the other hand, the Helfrich Hamiltonian is used to calculate the free energy cost due 
to a thermal fluctuation that changes the local curvature from its equilibrium value. The elastic coefficients are 
derived from the mean-square amplitudes of the fluctuations. The situation in which there exist two methods for 
calculating elastic moduli is reminiscent of other cases, for instance, the two different methods of evaluating the 
elastic constants of thermodynamic systems in linear elasticity theory pjjj, I40L El l42l | , and the two approaches 
for determining the surface tension of a planar interface |43t l44| . In the latter examples the different approaches 
lead to the same values for the mechanical moduli, in accord with the linear response i/ieori/|45ll46j. This is not 
the case with the rigidity constants of a liquid- vapor interface |37j . The discrepancy between the two methods 
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is due to the fact that in order to change the equilibrium radius of curvature of, say, a spherical liquid drop, it 
is necessary to change its volume as well. This means a change in the volume fractions of the two phases (i.e., 
the condensation of vapor or the evaporation of liquid) , and it thus requires the variation of the thermodynamic 
variables like the temperature or the chemical potential. In the fluctuation case the radius of curvature is varied 
by thermal fluctuations, while the thermodynamic variables are not altered. 

In this paper we discuss the statistical mechanics of fluid bilayer membranes. We derive expressions for the elastic 
coefficients a, k, and R of the membranes, relating them to the interactions and the correlation functions between 
the amphiphilcs forming the bilayer. We use these expressions for a Monte Carlo (MC) determination of the elastic 
coefficients of a bilayer membrane computer model. Unlike the expressions derived for the rigidity constants of a 
liquid-vapor interface, our expressions are such that they can be evaluated using a single MC run performed on the 
(quasi) flat membrane reference system only. This feature greatly simplifies the computational procedure, and makes it 
more efficient and well-controlled. Another important distinction between the membranes discussed in this paper and 
the system of liquid-vapor interface studied in Ref. j3f| is the fact that the mechanical and the fluctuation methods 
for determining their rigidity constants lead to similar results. Our expressions are derived using the mechanical 
approach, namely by calculating the free energy variations resulting from the change in the area and curvature of the 
membrane. The numerical values of the elastic coefficients which we obtain from these expressions are compared with 
the values extracted from a spectral analysis of the thermal fluctuations around the flat reference state. We find a 
very good agreement between the two methods. This agreement, which is expected by virtue of linear response theory 
(see discussion in section ^J, reflects the fact that the curvature of the membrane can be varied by changing the 
shape of the container (namely, by the application of external forces) without affecting the thermodynamic properties 
of the bulk aqueous phases surrounding it. It should be noted that the experimental values of k measured (for the 
same lipid bilayers) using mechanical and fluctuation methods can differ by as much as a factor of 3 |38l ] . The origin 
of these discrepancies is not well understood. 

The bilayer computer model which we use in this paper has been recently introduced by one of us |47j . (Here we 
use a slightly modified version of that model which we describe in section ItVl ) This model has two features which 
simplify the derivation of thermodynamic expressions for the elastic coefficients and the simulations performed for the 
calculation of these expressions. First, the simulations are conducted with no solvent present in the simulations cell, 
i.e., as if the membrane is in vacuum. This feature greatly reduces the number of atoms in the simulation cell, thus 
enabling us to simulate a relatively large membrane over a very long MC run. The ability to perform long MC runs is 
very important since the quantities whose thermal averages we try to evaluate are very "noisy" , and accurate results 
can be obtained only if they are measured for a large number of configurations. The other feature is the nature of 
the interactions between the molecules forming the membrane. In our computer model the amphiphilic molecules are 
modeled as trimers and the interactions between their constituent atoms are pair-wise additive. For such systems the 
derivation of expression for the elastic coefficients (see section II I H is simpler than for systems including many-body 
potentials. Our discussion in this paper is, therefore, restricted to central force systems only. 

The paper is organized in the following way: The theoretical aspects of our study are presented in sections [n] 
and IIIII In section [n] we describe the relation between the equilibrium and the fluctuation routes for determining 
the surface tension a and the bending modulus k of bilayer membranes, and explain why these methods (if used 
appropriately) lead to similar results. Then, in section ITTT1 we derive expressions for these quantities based on the 
equilibrium approach. Our expressions relate a and k to the interactions and the correlation functions between the 
"interaction sites" of the amphiphilic molecules. The numerical results are presented in section llVl where we calculate 
the elastic coefficients of our model system using the two methods and find a very good agreement between them. 
Some technical aspects of the simulations are discussed in the Appendix. Finally we conclude in section 



II. THE EQUILIBRIUM AND FLUCTUATION ROUTES TO MEMBRANE ELASTICITY 

Linear response is a fundamental theorem which relates the fluctuations of a system around its equilibrium state 
and the response of the system to weak perturbations . In the context of elasticity theory it provides a link 

between the shape fluctuations of thermodynamic systems and their elastic moduli. For example, when a 2D flat 
interface is slightly stretched or compressed from its equilibrium area Aq, the variation of the (small) surface pressure 
II is given by pLl!^ 

K A = -Mf A , (2) 



where A is the area of the interface and Ka is the stretching/compression modulus. The above relation provides one 
way to measure Ka- An alternative approach for measuring Ka is to consider the thermal fluctuations of the area A 
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around the equilibrium area Aq |40t |50( ■ The equipartition theorem suggests that in the low temperature limit when 
fluctuations around A are small 

{(A-A f)=^, (3) 

where ks is the Boltzmann constant and T is the temperature, while (• • • } denotes a thermal average. Linear response 
theory can be also applied to describe the normal, curvature-forming, fluctuations of the 2D interface. The discussion 
in this case (of normal fluctuations) is, however, somewhat more complicated. A proof of the equivalence between the 
equilibrium and the fluctuation routes to the surface tension a of a fluctuating interface had been presented with great 
clarity by Cai et al. |51j . Below we extend that proof and address the two routes to the bending modulus k as well. 
One important difference between the present discussion and the one presented in Ref. is related to the nature 
of the fluctuating surfaces in question. Here, we consider an elastic surface consisting of a fixed number of molecules 
whose area density is varied when it fluctuates. By contrast, the surface studied in Ref. [HlJ is incompressible and its 
area density is fixed to its equilibrium value. The variation of the total area of the latter is achieved via the exchange 
of molecules between the surface and the embedding solvent. A more detailed discussion on the differences between 
the elastic properties of compressible and incompressible surfaces appears in Ref. |Tfij . 

Let us consider a 2D surface that spans a planar frame of a total area A p which does not necessarily coincide with 
the equilibrium (Schulman) area Aq. The surface is free to undulate in the direction normal to frame. The ensemble 
of conformations which the surface attains is governed by a Hamiltonian TL (h(r)) relating the elastic energy to the 
conformation of the surface. The conformation of the surface is described by some "gauge" function h(f), where 
r = (x, y) label the points on the reference surface. The exact form of the Hamiltonian TL is unimportant and, in 
particular, it is not limited to the Helfrich form JJJ. As we are interested in moderately-fluctuating surfaces (with 
no overhangs), we shall use the the so called Monge gauge z = h(r), where h is the height of the surface above the 
frame reference plane. In what follows we will restrict our discussion to symmetric surfaces (such as bilayers) with no 
spontaneous curvature, i.e., with no preference to bend toward either the "upper" or "lower" side of the surface. In 
other words, we assume that the average conformation of the surface is flat and for each r 

{h(r))=0. (4) 

We also assume that the surface under consideration is mechanically stable, and that the validity of Eq. Q is not due 
to the partition of the configurations phase space into several sub spaces for which (h(r)) 0. 

If the frame (projected) area A p is not equal to the equilibrium area Aq then it is necessary to apply a tangential 
surface pressure in order to fix the area of the frame. If, in addition, normal forces are applied then relation breaks 
down. The function 

h(r) = (h(f)) (5) 

can be regarded as the strain field describing the deformed state of the surface. The free energy of a system subjected 
to a small deformation can be expanded in a power series in the strain variables. In full analogy to Hamiltonian l|T[l. 
we can write the Helfrich free energy of the surface in the following form: 

F(h)=F(h = 0)+o-(A (h) - A p ) + ^kJ 2 (h) + RK (h) + h.o.t, (6) 

where A (h) is total area of the surface defined by the function h(r), while J (h) and K (h) denote, respectively, the 
integrated total and Gaussian curvatures defined by 

P ee [ drJ 2 (h(r)) , (7) 

JA P 

and 

K = / drK (h(r)) . (8) 
Ja p 

In Eq.© we set the spontaneous curvature cq = [see Eqs.lQJ and iJSJ], and use the effective (normalized) values 
of the elastic coefficients which are different from the "bare" values appearing in the Hamiltonian (see discussion 
earlier in sectionP). The higher order terms (h.o.t) in Eq.© include both products of the small variables (A — A p )/A p , 
J 2 , and K, as well as terms involving the gradients of the local curvatures. The latter are assumed to be small since 
we consider only nearly-flat surfaces described by functions h which vary slowly in space. Since a, k, and R appear 
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as the coefficients of the free energy expansion in strain variables [Eq.©], they can be also related to the following 
partial derivatives 



dF 
II 



(9) 



h(r)=0 



d 2 F 

dJ 2 



(10) 



h(r)=0 



and 



dF 
dK 



(11) 



h(r)=0 



evaluated at the reference state h(r) = 0. 

Equations l j^) l- Hll(l express the equilibrium (mechanical) route to a, k, and R. The complementary fluctuations 
approach is more easily formulated in Fourier rather than in real space. Let us take a square frame (the reference 
surface) of linear size L p — y^Ap~, and discretized it into N 2 square cells ("patches") of linear size / = L p /N, where I is 
some microscopic length of the order of the size of the constituent molecules. Since the description of the membrane 
as a 2D continuous sheet breaks down on length scales below I, the surface has to be defined only over a discrete set 
of points {f g = (x g , y g )} each of which located in the center of a grid cell. Outside the frame region, the function can 
be defined by periodic extension of period L p , i.e. h (xg_ + n\L p , y g + n^Lp) = h (x g , y g ) where ri\ and ni are integer 
numbers. The Fourier transform of the (real) function h(r g ) is defined by 



hq - 



L, 



-iq-T g 



(12) 



where the two dimensional wave- vector q has N 2 discrete values satisfying 

{q x , q y = 2irm/L pi m = -N/2, N/2-1}. 

The inverse transform is given by 

I 



h(r g 



5Zv 



iq-r g 



(13) 



(14) 



If the topology of the surface is fixed and it does not form "handles" then the periodicity of the surface leads to the 
vanishing of the Gaussian curvature (JSJ (Gauss-Bonnet theorem). Writing the expressions for the area A (h/j and the 
integrated total curvature J in terms of Fourier coordinates: 



A(h)=A p + -J2 q 2 hfh^+0 (|V 



(15) 



and 



J 2 (h) = I 2 tfhh-q + O {\h 9 \ 



(16) 



and substituting them in Eq.l|S|l. we obtain the following expression for the free energy 



F(h)=F{h = 0) + -J2 W + «9 4 



0(q 6 )] hfh. ? +0(\h q ^). 



(17) 



The free energy Ijl7(l can be related to the surface Hamiltonian 7i ({h(f g )}) via the partition function. We may use 
the Fourier transform 



hq- 



(18) 
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of the function h ({f g }), and express the Hamiltonian as a function of the Fourier modes: H ({hq-}). Introducing the 
set of Lagrange multipliers {jq-} each of which enforcing the value of hq- — (hq-) , we write the partition function of the 
surface as 



Zg [A p , {if}] = J D [{hf}] exp J -[3 



H({j?})-J2 h tft 



where j3 = (k B T) 1 . The associated Gibbs free energy is 

G[A p ,{j 9 }] = -k B T\nZ G . 
From Eas. l|19|) and l|2(J|) it is easy to derive the following relation 



hq = (hq) 



dG 



and 



d 2 G 

{hg-h-q) - (/lf)(/l-f) = -k B T -jT-jj - . 

The Helmholtz free energy F is related to G via 

F [A p , {%}] = G [A p , 0»] + hk> 



(19) 



(20) 



(21) 



(22) 



(23) 



where 



dF 
dhg- 



Jq- 



If we use expression <|17[) for the Helmholtz free energy, we find from Ea. l24fl that 

i f = I 2 [aq 2 + nq 4 + 0(q 6 )] +■■■ 



(24) 



(25) 



[note that hq- (jq* = 0) = 0]. Combining Eas. lfTTjl . l|2*3*|l. and (|2*5|l we obtain to the following expression for Gibbs free 
energy 



G = F({h q j={Q})-Y; 



Jq.J — q 



2l 2 [aq 2 + Kq 4 + 0(q e )} 



(26) 



When this expression for G is substituted in Ea. (|22|) and evaluated for {jq-} = {0} (which corresponds to the reference 
state {hq*} — {0}), we find that the mean square amplitude of the fluctuations with a wave- vector q (the "spectral 
intensity") is given by 



(hq*h-q*) 



{M={o} 



\hq\ 2 ) 



k B T 



{M={0} 



nq* 



(27) 



This result, which quantifies the magnitude of the fluctuations around the flat equilibrium state, provides a second 
( "fluctuation" ) route for calculating a and n (but not for the saddle-splay modulus R) . It is frequently quoted in an 
incorrect form with do and ko, the coefficients in the Helfrich Hamiltonian instead of a and k. The equivalence 
of the two routes to membrane elasticity is expressed by the fact that the elastic coefficients appearing in expression 
(|27H are the same as those obtained from Eas.(|^l— l|lll) . and which are associated the "equilibrium" route. In the next 
section we use Eas.l|9T)- (lllfl to derive statistical-mechanical expressions for the elastic coefficients. Then, in section llVl 
we demonstrate, using computer simulations of a bilayer membrane model, the agreement between the two different 
methods of calculation. 
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III. THERMODYNAMIC EXPRESSIONS FOR THE ELASTIC COEFFICIENTS 

A. The surface tension 

Let us return to the equilibrium route to membrane elasticity and to expressions ijHJ which describe the relation 
between the free energy and the elastic coefficients. The surface tension can be computed by comparing the free energy 
of the membrane at the reference state (which is assumed to be flat) and the free energy of a flat membrane with 
a slightly larger area. These two membranes are shown schematically without the underlying microscopic details, 
in Figs, n (a) aud (b). We reemphasize that the total number of amphiphilic molecules which form the membrane 
is fixed, and that the surface tension should be related to the free energy dependence on the area density of the 
amphiphiles (rather than the free energy cost to add molecules to the membrane). The characteristic surface of the 
membrane to which the free energy is applied, is chosen as the mid surface between the two layers. The total volume 
of the membrane is assumed to be fixed; otherwise, an additional term involving the volume compression modulus 
must be introduced in Eq.©. 

It is important to remember that in Figs.^( a ) an d (b), only the mean configurations of the surface (in the reference 
and deformed states) are depicted, and that the surface undulates around these (ensemble) average conformations. 
In other words, "the state of the surface" refers to its average conformation. As has been discussed earlier in section 
ITT1 normal forces must be applied in order to deform the surface from its reference state |52j . If the membrane is 
embedded in a solution and placed in a container, than these forces can be generated by deforming the entire container, 
as demonstrated in Fig.^(d). Such a system can be conceptually divided into bulk aqueous phases and the interface 
between them which includes the membrane and the adjacent hydration layers. The volumes of the bulk phases above 
and below the membrane are fixed by the presence of solute particles that cannot permeate the membrane. The 
deformation of the boundaries of the container "percolates" to the interface and the latter acquires the shape of the 
surface of the container. However, since the bulk solution is fluid and has a vanishing shear modulus, its deformation 
without changing its volume does not add any contribution to the free energy. 

Even thought real bilayer systems are always embedded in a solvent (which influences their elastic properties), the 
calculation of the surface tension can be also performed for model systems that exclude the latter and leave only 
the interfacial region. This is possible due to the fact that the surface tension can be calculated by considering a 
deformed flat membrane. Such a membrane can be uniquely defined by the perimeter V (h(r)) of the characteristic 
surface [represented by open circles in Fig.^(d)]. The free energy of the membrane can be derived from the partition 
function Z via the relation 

F = -k B T\nZ. (28) 

The expression for the partition function must take into account the microscopic nature of the membrane, and the 
potential energy E due to the interactions between the amphiphilic molecules. In what follows we assume that E can 
be written as the sum of pair interactions between the atoms ( "interaction sites" ) forming the molecules 

E=Y,<t>{r^), (29) 

(a/3) 

where the distance between atoms a and /3, and summation over all pairs of atoms {a(3) is performed. The 

various interactions are not identical but rather pair-dependent, as each amphiphilic molecule is typically composed 
of many different atoms. They also depend on whether the atoms belong to different molecules or part of the same 
amphiphile. In the latter case some atoms are covalently bonded what brings in an additional contribution to E. For 
brevity we will omit the subscripts of the potential and the indices of the argument r a P will serve as an indicator of 
the specific potential. With the potential energy described by Eq.J^HJl, the partition function is given by 

Z = Yl «p(-5> (r^) /k B T I , (30) 

V - Conf. \ (af3) / 

where the sum runs over all the conformations is which the perimeter of the characteristic surface is depicted by 
the closed curve V '. Our assumption that the membrane has no spontaneous curvature guarantees that its average 
conformation is indeed flat. Alternatively, one may consider the system together with the bulk phases, and replace 
the sum in Ea. H3U|) with integration of the coordinates of all atoms {r 7 } over the entire volume of container (or the 
simulation cell) V^ n 
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In addition to the above integral, it is necessary to specify the boundary conditions for the positions of the amphiphiles 
near the walls of the container, so that the perimeter of the characteristic surface would be described by V. 

Let us assume that our cell (container) has a square cross section of linear size L p with —L. p /2 < x,y < +L p /2. 
The deformation of the cell depicted in Fig. ^ (d) can be described by the following linear transformation 





(32) 



which maps every point R on the boundaries of the undeformed cell to its strained spatial position r. The characteristic 
surface has the same shape as the upper and lower faces of the cell, and its area is given by 



A = A P VT 



An 1 



(33) 



where A v = L p is the area of the reference surface. Since the deformed surface which we consider is flat, its free 



energy is given by [see Eqs.© and (|3*5)l ] 



F = F(e = 0)+aA p - 



(34) 



from which we conclude that 



a = 



1 cPF 
A p de 2 



(35) 



e=0 



Using the relation between the free energy and the partition function (|28|) , we may also write the above result in the 
following form 



1 <Pz 
z'dT 2 ' 



i 

^2 



dZY 
Ik) 



(36) 



e=0 



If we now turn to our expression (|31|l for the partition function, we notice that it depends on e only through 
the integration volume V^ G n- The differentiation of Z with respect to e, however, could be carried out more easily 
if the dependence on e is removed from V cc u and brought into the integrand. In other words, we wish to change 
the integration variables in l|31(l from r 7 to R 1 , where the latter are confined inside the undeformed cell. This is 
achieved using transformation (|32|) . which originally described the deformation of the boundary points, and is now 
being applied inside the volume of integration [53j. With the new set of variables, the distance between two atoms is 
given by 



r «/3 



R af} ) +2eRfR 



a/3 pa/3 



e 2 (i? 



a(3\ 



1/2 



(37) 



In the undeformed reference state r a/3 (e = 0) = R a<3 . The partition function reads 



Z 



JV ° 7 = 1 V < Q ' 3 ) 



(R a Pf + 2eRfR« P + e- 



(Rf)' 



/k B T 



(38) 



where Vq = V ce \i(e = 0) is the volume of the undeformed cell. The Jacobian of the transformation has been eliminated 
from the integrand in the above expression since it is unity. The differentiation of Z with respect to e is now 
straightforward but lengthy. We skip the details of the calculation, and write below the final expressions for the first 
and second derivatives, evaluated for e = [only the value at e = is required in Ea. (|36(l ] 
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and 
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where (/>' = d(f>/dR and 0" = d 2 4>/dR 2 . When these expression are substituted into Ea. (|36fl we readily find that 
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where the thermal averages are evaluated at the undeformed reference state of the system (e = 0). If the system 
is macroscopically invariant with respect to reversal of the sign of the z coordinates (z — > —z), then the first term 
in the above expression for a vanishes. If, in addition, the system is invariant with respect to rotation around the 
z axis (x — > y; y — + —x), then another expression can be derived with R"@ replaced by R"E ! . Defining R 



a 6 



(jlx^j + (Ry^ j we finally arrive to the following expression: 
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This expression can be also written in the following compact form 

a = L 



CxZXZ H~ Cy Z y Z PxX 



L z /d z 



(42) 



(43) 



where L z is the linear size of the system (the cell) in the z direction (normal to the membrane), while P and C denote 
the pressure tensor and the tensor of elastic constants of the system. The quantity fi z t is the shear modulus associated 
with the deformation depicted at Fig.Q](d) 0, Ell- 
in is interesting to compare the above results H42l) - (|43[) with the much better known (and more frequently used) 
expression for the surface tension l34| 

1 M'-2(Rf)\_ 
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(44) 



The latter expression is obtained when one considers the variation of the free energy resulting from the (volume- 
preserving) variation of the projected area A p 
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(45) 
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The deformed state associated with the surface tension a is shown in Fig. Q](c). For fluid membranes we anticipate 
that a — a since the difference between them 



C x 



a 



yzyz 



(46) 



is proportional to the shear modulus fXtz associated with the deformation shown in Fig. ^ (e) . The shear modulus 
fitz is expected to vanish because the areas of the characteristic surfaces of the membranes in Figs. ^ ( a ) an d (e) are 
identical; and the Helfrich free energy of a flat membrane depends only on the area of the characteristic surface, but 
not on the orientation of the plane of the membrane with respect to the walls of the container. This argument for 
the coincidence of a and a could be applied directly to the membranes in Figs. ^ (b) and (c), whose characteristic 
areas (as well as their volumes) are also identical. The tilt of the cell's wall, however, can be safely ignored only in 
the thermodynamic limit, when the width of the membrane becomes much smaller than its lateral dimensions. If the 
system is not sufficiently large than the Helfrich form for the free energy in which the membrane is associated with a 
2D characteristic surface is not entirely applicable. The finite width of the membrane must show up in the expression 
for the free energy, and the surface tensions a and a do not perfectly agree. 



B. The bending modulus 

The bending modulus can be calculated by considering a deformation of the characteristic surface from a flat to 
cylindrical geometry. The deformation, which is depicted in Fig. [21 can be described by the following nonlinear 
transformation of the boundaries of the cell [compare with Eq.©] 



T x — R x 

Ty = Ry 

r z = R z 



R l 



Rl Lj/4, 



(47) 



where —L p /2 < x < +L p /2, and i?o 3> L p is the radius of curvature of the cylinder. The integrated total curvature 
[Eq.Q] of the characteristic surface is 



J = — 

Rq 



(48) 



and its area is 



A = A„ + 2 arcsin 
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(49) 



The free energy is, hence, given by 
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(50) 



from which we deduce the following relation 
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where 



Rq 



(52) 



The calculation of the r.h.s of the above equation is very similar to the one presented in section ITlI Al which was 
based on expression (|38(l for the partition function. The deformed pair distance, which in that case was given by 
Eci.l|H7jl. is now depicted by the following relation 
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(53) 
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where 



X af} 



x a + x fj 



(54) 



is the average of the x-coordinates of atoms a and /3 (in the undeformed state). Comparing Eas. i|37|) and i|53|l . and 
respectively, Eas. (|35|) and (|51|l . it is easy to realize that the result of the calculation is the following expression 
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(55) 



\( a {3) 



which is similar to Ea. H41|) . except for the fact that R^P is everywhere replaced by —X a ^R^P. 

Among the five terms on the r.h.s of Eq.j55j), only the second involves averages of quantities including the product 
X a/3 X lS with (aj3) ^ {jS). In the other four terms, the quantities X a P and (X Q/3 ) can be replaced by their averages 
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(56) 



(57) 



since they multiply quantities which depend only on the separation between atoms a and (3 and whose averages, 
therefore, arc independent of the location of the pair (provided the system is invariant to translations in the x and y 
directions). This, in combination with Ea. H41|) . yield the following expression for n 
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(58) 



Replacing R%P with R^ , and X a0 with Y a(3 = (Y a + Y )/2, we obtain the "symmetric" formula 
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(59) 



It is important to remember here that the above expression for k (I59|l applies to square membranes only with the 
origin of axes located at the center of the membrane so that —L p /2 < x,y < +L p /2. A formula which does not 
depend neither on the shape of the membrane nor on the location of the origin is obtained as follows: The first and 
third terms in Ea. (|59J) can be written jointly in the following form 
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and 



where 

jpa/9, 7 * = ^ (62) 

va.8 viS 

ee * ~* ■ (63) 

The terms appearing before the square brackets in Ea. H61|) depend only on the relative coordinates of atoms with 
respect to each other. Therefore, the average of (i a,S ' li ) (the second term in square brackets, which depends only 
the location of the center of the pair/triplet/quartet in question) can be performed separately. As in Ea. H57|) we 

have ^(X a/3,7<5 ) 2 ^ = L^/12, what leads to the cancellation of the first two terms in square brackets in Eci. (|61[) . 
Applying the same argument for the second and fourth terms in Eq. I|59|) , and defining 

F a/3, 7 «5 = (Y a P + Y~< 5 ) /2, and 

Ay' 3,7 ' 5 = (Y a @ — Y l5 ^ /2, we arrive to the following expression 



2A„k B T 



(64) 



which is the more general form for expression (|59|l since it is independent of the shape of the membrane and of the 
location of the origin of axes. 

The deformed membrane, shown schematically in gray shade in Fig. [5J may be considered as part of a closed 
cylindrical vesicle (depicted by the dashed line Fig. EJ) . Accordingly, one may argue that its free energy is given by 

6 

F — — inside (65) 

where -F V csicic is the free energy of the vesicle and 9 is the apex angle of the deformed membrane. This relation, 
however, is incorrect since -F V esicic includes a term which is unique to closed vesicles and should be omitted in the 
case of open membranes. The additional contribution to fVesiele which has been termed "the area- difference elastic 
energy" , should not be confused with the bending energy. The latter is the free energy required to bend the membrane 
while keeping its area density fixed. The former, on the other hand, originates from the simple fact that upon closure 
of the vesicle, it becomes impossible to preserve the area densities of the amphiphiles in both the outer and the inner 
monolayers. The outer monolayer is stretched and the inner monolayer is compressed relative to the mid characteristic 
surface. The elastic energy resulting from such curvature-induced changes in the monolayer areas is a non-local effect 
because the monolayers are capable of independent lateral redistribution to equalize the area per molecule of each 
leaflet. The distinction between (local) bending elasticity and (non-local) area-difference elasticity has been discussed 
by Helfrich, not long after introducing his famous Hamiltonian j54| . The idea, however, did not gain much popularity 
until the issue has been analyzed systematically by Svetina et al. some years later j55j. Early theoretical works and 
experimental measurements of the bending modulus failed to separate the local and non-local contributions |56j . This 
is not the case with our expression l|64|l for k which has been derived by considering an open membrane. For an 
open membrane, the two leaflets have the same area as the top (button) surface of the containers and, consequently, 
area-difference elasticity do not show up. 



C. The saddle-splay modulus 

Finally, we derive our expression for the saddle-splay modulus R. The following transformation, applied to the 
boundaries of the container 



r z = R z + yjl% - - y 2 - yj R 2 - (66) 

(with —Lp/2 < x < +L p /2), describes a deformation of the surface to spherical geometry where the sphere's radius 
Rq 3> L p . It is not difficult to show that the free energy of the spherical surface is given by 

H 2 + --- , (67) 
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where H = 1/Rq. From the above expression for F, the following relation 
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is easily derived. The deformed pair distance is 
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where X a ^ and Y a/3 have been defined in section ITlIBI Since Eas. (|68|l and (|69|l have, respectively, the same form as 
Eqs.JSlJ and (|53(l . we immediately conclude that the r.h.s of Ea. (|68|l is given by expression similar to (|55|> in which 
X a PR<f is everywhere exchanged with X^R'f + Y af3 R^. Following the same steps described in the derivation of 
Ea. (|59|l from l|55|) . and using the additional relation 



we finally arrive to the following result 



L p /2 r L p /2 



T 2 , 

J ~L p /2 J -L p /2 



ApksT 



5>' (R a " 

<a/3> 



X aSi R<fR^ 
RPV 



xy dxdy = 0, 



yafi r>af3 r>a/3 

(a/3) 



(70) 



(71) 



This expression applies to square membranes only, with the origin located at the center of the membrane. The more 
general expression is 



K = —K — 
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(72) 



IV. NUMERICAL RESULTS 



The purpose of the MC simulations which we conducted and present in this section is twofold: The first is to 
test the validity and accuracy of our expressions for the elastic coefficients. The second is to examine the agreement 
between the mechanical and the fluctuation routes to membrane elasticity, as discussed in section [H] The model 
system whose elastic properties were studied by the simulations has been described in great details in Ref. |47j . 
Briefly, the "lipids" that serve as the building blocks of the membrane consist of three spherical atoms of diameter 
a (see Fig. 0) interacting with each other via pair- wise L J potentials (whose details can be found in Ref. ) ■ To 
avoid the complications involved with long-range interactions, the LJ potentials have been truncated at some cut-off 
separation R a/3 = r c = 2.5a and, in addition, modified to ensure the vanishing of <j> and its first two derivatives, 
0' and <j)" , at r c . The continuity of the second derivative of the pair potentials is an important feature since <fi" 
appears in our expressions l|42|l for a . Two changes have been made in comparison to the original model presented 
in Ref. |47l| . The first is a small reduction of the temperature which, in this paper, has been set to 0.9To where To is 
the original temperature (in Ref. |47^1. The second is the addition of new interactions between atoms which are part 
of the same molecule. In Ref. |47| the molecule were linear rigid trimers with a fixed distance a between the centers 
of the constituent atoms. Here, we allow some little variations of the separation between the atoms. The mid atom 
(labeled 2) has been linked to the two end atoms (labeled 1 and 3) via harmonic springs with spring constant K and 
equilibrium length a: 

0(R) = ^K(R-a) 2 , (73) 

while the pair potential between the end atoms has been set to 

<t>(R) = l -K{R - 2a) 2 . (74) 

We use a large value for the spring constant K = 8000 kgT/ a 2 , for which the separations between the atoms do not 
exceed the order of 1% of their equlibrium values. While this means that the molecules in our model are "almost" 
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linear and rigid, the use of the above potentials (|73|l and (|74(l creates a situation in which all inter-atomic interactions 
(whether between atoms belonging to the same or different molecules) are depicted by smooth potentials; and so, 
our expressions for the elastic constants can be used without any further complications. The total number of lipids 
in our simulations was N — 1000 (500 lipids in each monolayer), and no additional solvent molecules were included 
inside the simulation cell (as if the membrane is vacuum). Periodic boundary conditions were applied in the plane of 
the membrane, and no boundaries for the simulation cell were defined in the normal direction. The linear size of the 
(square) membrane was set to L p = 29.375a. Subsequent MC configurations were generated by two types of move 
attempts: translations of lipids (which also included some minute changes in the relative locations of the three atoms 
with respect to each other) and rotations around the mid atom. A set of 2N = 2000 move attempts of randomly 
chosen molecules is defined as the MC time unit. Both types of moves (translations and rotations) were attempted 
with equal probability, and the acceptance probabilities of both of them was approximately half. The MC relaxation 
time has been evaluated in Ref. ■ It is of the order of 10 4 MC time units and has been very little affected by the 
changes introduced in the model. A typical equilibrium configuration of the membrane is shown in Fig. 



A. The fluctuation route 



The fluctuation approach for determining the surface tension a and the bending modulus k is straightforward to 
implement: The profile of the membrane in our simulations was defined by mapping the system onto an 8 x 8 grid, and 
defining the height h({f g }) of the membrane in each grid cell as the average of the local heights of the two monolayers. 
The latter were evaluated by the mean height of the lipids (whose positions were identified with the coordinates of 
their mid atoms) belonging to each layer, which were instantaneously located inside the local grid cell. Note that the 
mesh size I — L p /8 ~ 3.67a is somewhat larger than the size of the lipids, as required in our discussion in section 
ITT1 The Fourier transform of h({r g }) was obtained using Ea. (|18|l . and the mean squared amplitudes of the different 
modes were, eventually, fitted to the inverse form of Ea.i(2"7)l 

1 = l 2 [aq 2 + K q 4 + Q(q 6 )} 
(N 2 ) k B T [ > 

The results of this spectral analysis are summarized in Fig. [5] where we plot the value of l/l 2 (\h^\ 2 ) as a function 
of q 2 . The error bars represent one standard deviation in the estimates of the averages, which were obtained from 
simulations of 16 different membranes and a total number of 1.25 x 10 4 measurements of the spectrum per membrane. 
The measurements were done at time intervals of 100 MC time units. The curve depicts the best fit to Eq.J7SJ), which 
is obtained when a and k take the following values: 

, knT 

a = (-0.6 ±0.2 -V 

a z 

k = (46 ±2) k B T. (76) 

The contribution of the g 6 term to the fit was, indeed, significantly smaller than that of the other two terms on the 
r.h.s. of Ea. lf75|l . 



B. The equilibrium route 



While the measurement of a and k using the fluctuation approach was a relatively straightforward matter, the 
application of the equilibrium approach emerged as somewhat more challenging task. The most significant differences 
between the two approaches was the amount of computer resources required for an accurate determination of the 
elastic coefficients. The results which we present in this section have been obtained using 64 nodes on a Beowolf 
cluster consisting of Intel architecture PCs, where the CPU time per node was of the order of three months. The need 
of such a large computer time should be compared to the relative ease with which the results in Eq. I|76|l have been 
obtained - using a total number of only 16 nodes over a period of about 10 days. The reason that the equilibrium 
approach is so much computer-time-consuming is the "noisy" nature of the statistics of the terms whose averages are 
evaluated in expressions (|42|l and (|59|l . From the conceptual point of view, the determination of the surface tension 
a using expression (|42|l is pretty simple. The determination of the surface tension a from expression (|44l) is even 
easier since it is a much less noisy quantity. In fact, the computational effort required for an accurate determination 
of the value of a is even smaller than the one required for the calculation of a by the fluctuation method. The surface 
tension a does not apply directly to membranes with a fixed projected area. Yet, it is expected to coincide with a in 
the thermodynamic limit 



15 



The determination of k is more complicated. Here we can, in principle, choose between expressions (I59|l and i|64[l . 
The latter is more general (since it is not restricted to square membranes), but prohibitively time consuming. This 
can be understood by considering the number of operations required for a single measurement of the quantities of 
interest. Assuming each atom in our simulations interact with a finite number of other atoms, the total number of 
operations required by expression l|64l) is 0{N 2 ), while the number required by expression l|59[) is only O(N). In our 
simulations the total number of atoms is 3000, which means a difference of about 4 orders of magnitude in efficiency 
Using expression (|59|l to measure k is, however, tricky because this expression involves not only the relative locations 
of the particles with respect to each other (as in the case of the expressions for the surface tension), but also the 
absolute coordinates of atoms. This would not create a problem if only the central coordinates (X a P and Y a @) of 
the pairs had to be found [as one may, naively, conclude from Ea. l|5^|l ]. since that among the set including the pair 
(a, (3) and all its periodic images, only one satisfies the requirement —L p /2 < X a>3 , Y a @ < +L p /2. However [and this 
becomes clear from the derivation of expression i|64fl from expression l|59f) ] , what we actually have here is a periodic 
boundary conditions problem where the pairs play the role of the particles, and X a @ and Y a P serve as the coordinates 
of these "particles" . This means that each quartet ((a, /3), (7, S)) is identified as the pair (a, (3) and the pair (7, 5) or its 
image nearest to (a,/3) and, in addition, that the center of the quartet must satisfy —L p /2 < A >a/3,7<5 , Y a ^' jS < L p /2. 
The fact that sometimes a pair must be replaced by one of its images (which are located outside the boundaries of 
the simulation cell) is problematic since this means that the location of the pair, which is needed in expression i|59[l . 
cannot be specified by a single value. A solution to this problem is obtained by dividing the simulation cell into 
stripes parallel to either the x or the y axes [depending on whether we calculate the third or fourth term in Ea. H59|) ]. 
and to split the summation over all the pairs to several partial sums over the pairs included in the different stripes. 
The partial sums corresponding to the images of each stripe (which consist of all the images of the pairs included in 
the stripe) can be found with almost no additional effort. The product of two partial sums gives the contribution of 
all the quartets consisting of pairs located inside the two relevant stripes. Depending on the distance between the 
stripes (along the relevant axis) and their locations with respect to the center of the cell, it is usually easy to decide 
in which case a stripe should be replaced by one of its images. Ambiguities about the correct decision occur in a finite 
number of cases (i.e., for a finite number of pairs of stripes). In these cases, individual decisions must be made for 
each quartet. The number of such quartets can be reduced significantly if the system is divided into a large number of 
stripes N s , since the narrower the stripes the smaller the number of pairs included in each one of them. A more elegant 
solution is to choose a certain convention about the ways the contribution from the ambiguous quartets is added to 
Eq. (|59|l . This will inevitably introduce a systematic error to our estimates of n. However, if we make a set of estimates 
based on increasingly larger values of N s , we can obtain the correct averages by extrapolating our results to the limit 
1/A^s — > 0. The method, which is described in more details in the Appendix, can be generalized to handle correctly 
the calculation of R. However, because of the mixing of the x and y coordinates in Ea. H71l) . the implementation of 
the method becomes more complicated. For this reason, and due to the fact that the fluctuation approach does not 
provide a value of saddle-splay modulus to compare with, we did not use our simulations to calculate R. 

In section[H]we have explained in great details why the elastic coefficients obtained from the fluctuation approach are 
the free energy coefficients a and k rather than the Hamiltonian coefficients Co and k,q. This means that the quantities 
in expressions l|42|) and (|59|l should be averaged over the ensemble of all possible microscopic configurations. However, 
it is also easy to understand that the same expressions can be used to calculate the Hamiltonian coefficients. The 
latter, which characterize the energy changes caused by deformations of the flat membrane, can be obtained by 
restricting the averages to conformations where h(f g ) = for every grid cell, thus avoiding the entropic contribution 
of the thermal fluctuation to the free energy. To sample this configuration phase space one need to accompany every 
MC move attempt with one or two (depending on whether the molecule leaves the grid cell or not) additional moves 
of adjacent molecules. Moreover, one can also sample the phase-space consisting of only those conformations of the 
membrane with wave vectors in the range 2ir/L p < A < q. The results of such a calculation are the wave-dependent 
coefficients cr(A) and k(A). One of the problems which can be studied by such inve stig ation is the value of the 
numerical coefficient c in the formula for the renormalized bending modulus [53, IH, |5|| |6(| : 

k(A) = k + c^MAI). (77) 

47T 

This problem aroused a renewed interest recently since its has been suggested that the value of c may be positive, 
which means (quite remarkably) that the fluctuations stiffen rather then soften the membrane |6ll 163 . 631. 

While determining the value of c was not possible with our computer resources, we did use Ea. (|77|) in our analysis 
of the results. Our need of Eci. (|77Jl and the link that it provides between k and ko is related to the peculiar nature of 
our simulations which are made in a "solvent-free" environment. As has been discussed in section lHll our expressions 
for the elastic coefficients have been derived based on the assumption that the membrane is embedded in solvent and 
that the entire container is deformed. In our simulations, however, we have no container (there are no boundaries 
for the simulation cell in the z direction) and, so, the applicability of our approach should be examined carefully. 
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The arguments which we presented in section IlII Al [see, in particular, the discussion around Ea. i|30[l ] demonstrate 
that the presence of solvent is essential only for the calculation of k and R, but not for the calculation of the surface 
tensions a and a. By contrast, the Hamiltonian coefficients can be all measured in a "solvent-free" model since they 
are extracted from simulations of flat, non-fluctuating, membranes. The value of kq and the relation given by Eq.J77|) 
provide then an estimate for the value of k. Since the finite-size correction to the value of k grows only logarithmically 
with the size of the system, and since kq S> fcsT, the difference between k,q and re is not significantly large. In our 
simulations it actually falls within the uncertainty in our estimates of the bending modulus, which means that k and 
kq are practically indistinguishable. In addition to our measurement of Ko, we also measured k directly from the 
simulations. As we have just explained above, such a measurement is expected to fail and to lead to the incorrect 
conclusion that k = 0. We used this incorrect result as a test for our code. 

The values of the elastic coefficients have been extracted from simulations of 64 membranes starting at different 
initial configurations. The initial configurations were generated by randomly placing 500 lipids in two different layers 
with a vertical (along the z direction) separation a (the size of the atoms) between them. The initials configurations 
were "thermalized" over a period of 2 x 10 5 MC time units, followed by a longer period of 1.2 x 10 6 time units 
during which quantities of interest were evaluated. The uncertainties in our final results correspond to one standard 
deviation in the estimates of the averages. We first made the simulations with non-fluctuating membranes, from 
which we extracted the values of the Hamiltonian coefficients. Then, we removed the part in our algorithm which 
is responsible for keeping the membrane flat. The membranes were equilibrated again, and then the values of the 
thermodynamic (free energy) coefficients were determined. 

For the bare coefficients we find the following values for the surface tension: 

a Q = (0.8 ±0.5) 

ao = (-0.07 ±0.01) (78) 
a 1 

The comparison of these results with each other, and with the values of the elastic coefficients extracted from the 
fluctuation approach [Eq. l|76|l ] reveals: (a) a disagreement between the two surface tensions ao and ao, which should 
be attributed to the finite size of our membrane (see our discussion in section Till A(l : and (b) a disagreement between 
ao and a which should be attributed to the entropic contribution to the surface tension. The bending modulus kq 
has been obtained by dividing the system into N s stripes and extrapolating the results for ko to the limit 1/N S — > 0, 
as explained earlier in this section (see also the Appendix). From the extrapolation procedure, which is summarized 
in Fig. El we find that 

n = (44 ±10) k B T. (79) 

This result also serves as our estimate for k (see discussion earlier in this section). The similarity of the above value 
of k (which is, unfortunately, obtained with a somewhat large numerical uncertainty) to the one quoted in Ea. H76|) 
corroborates the argument presented in section [n] regarding the equivalence of the two routes to membrane elasticity. 
Further support to this argument is obtained from the agreement of our result in Ea. l|76|l to a, with the value of the 
surface tension obtained from equilibrium approach [using expression (|42|) ]: 

a= (-0.3 ±0.5) (80) 
er 

Our result for a [expression l|44|l ] is not very much different 

a= (-0.41 ±0.01) (81) 

These values are quite different from those given in Eq. l|78|l . thus demonstrating the importance of the entropic 
contribution to the surface tension. 

Finally, we plot in Fig. [7] our results for the "apparent" bending modulus k* which we have obtained, using 
expression l|59[) . from simulations of a fluctuating membrane. These simulations serve as a test for our code. We find 
k* = (—4 ± 8) ksT which is consistent with the anticipated value n* = 0. 



V. SUMMARY AND DISCUSSION 



Motivated by the lack of a well accepted theory to deal with the statistical-mechanical behavior of curved interfaces, 
we have studied the elastic properties of fluid bilayer membranes using analytical and computational tools. Two 
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distinct methods have been employed to measure the surface tension a, and the bending modulus k, of a model 
membrane. In the first ("fluctuation") method the elastic coefficients were extracted from the analysis of the spectrum 
of thermal fluctuations of the membrane. The second ("equilibrium") method is based on the fact that a and k describe 
the free energy variations due to area-changing and curvature-forming deformations and, therefore, can be related to 
the derivatives of the partition function with respect to the relevant strain variables. Using this kind of relation, we 
have derived formal expressions for a and k in central force systems. Our expressions associate the elastic coefficients 
to the interactions between the molecules and the two-, three-, and four-particles distribution functions. The most 
important feature of these expressions is the fact that even though a and k (as well as the saddle-splay modulus R) 
are related to deformations of the membrane, they can be extracted from a single MC run performed on the reference 
(unstrained) system. 

One of the puzzles about curved interfaces elasticity is related to the correspondence between the above two 
approaches for determining their rigidity constants. We used linear response theory to prove that the two methods 
must agree for the values of a and k provided that the system is deformed by the application of external forces and 
not by altering other thermodynamic variables such as the temperature or the chemical potential of surface molecules. 
Moreover, our discussion clarifies that the coefficients in question, a and k, are the effective elastic coefficients which 
appear in the Helfrich free energy (rather than the Helfrich Hamiltonian) and which are influenced by the thermal 
undulations of the membrane. Our computer simulations and the numerical values of the elastic coefficients which we 
find, confirm the idea of equivalence between the two routes to membrane elasticity. 

Comparison of the computational efficiency of the two methods shows that for our membrane model system the 
fluctuation method provides more accurate estimates of the elastic coefficients than the equilibrium method, and 
requires less CPU time. The major shortcomings of the fluctuation approach is the fact that it can be utilized for 
measurements of the effective coefficients only, and that it requires the determination of the profile of the interface 
during the course of the simulations. While this is easy with our "water-free" computer model, this may not be so 
in other cases, for instance, for membranes which tend to exchange molecules with the embedding solvent, or for 
liquid-vapor interfaces near the critical point when the interface is difficult to distinguish from the bulk phases. In 
these cases the equilibrium method may be more attractive since the interactions in the bulk phases do not contribute 
to the values of a and k when calculated using our expressions for the elastic coefficients. Moreover, with the same 
mechanical expressions for a and k, the bare (Hamiltonian) coefficients can be also calculate. Our measurements 
demonstrate that close to the tensionless state of the membrane, the entropic component of the surface tension is 
quite significant. This has been also found recently in a theoretical study of the surface tension of fluctuating surfaces 

Finally, we would like to reemphasize that our expressions for the elastic coefficients apply for central force systems 
only. Following our derivation of these expression one should be able to generalized them to more complicated cases 
including many-body interactions. A more realistic model must also include electrostatic interactions whose long-range 
nature pose a computational challenge. 
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APPENDIX A: DETERMINATION OF k USING THE METHOD OF STRIPES 

The most common way to reduce finite size effects in computer simulations is obtained by employing periodic 
boundary conditions, namely by regarding the simulation cell as part of an infinite periodic lattice of identical cells. 
When the range of the interactions is less than L p /2 (half the linear size of the cell) than each particle a interacts only 
with the nearest periodic image of any other particle 0. This, in turn, is identified as the pair (a,/3). Each pair has 
infinitely many periodic images each of which is associated with a different simulation cell; and with each simulation 
cell each pair is associated exactly once. The set of all the different pairs associated with one of the cells [say, the 
original ( "primitive" ) cell] is the one over which the summation in expressions l|42|l and l|44|) for the surface tension 
should be performed. 

Things become more complicated when we try to evaluate the bending modulus k using expression l|59|) . In this 
case, coordinates associated with the location of the pair {X a P and Y a P) appear in the expression, and so it becomes 
necessary to decide which of the periodic images of each pair is actually associated with primitive simulation cell 
{—L p /2 < x,y < +L p /2) over which the sum in Ea. (|59|l is performed. The intuitive candidate is the periodic 
image with —L p /2 < X a @ ,Y al3 < +L p /2. Making this choice, however, is not the right convention. The correct 
way to handle the summation in expression (|59|l can be deduced from our derivation of expression i|tj4[l which is 
independent of the location of the origin of axes. Following the discussion that led from Eq.JS^J) to Ea. (|64|) it becomes 
clear that: (a) each quartet ((a, (3), (7, 5)) must be reproduced exactly twice from sums in Ea. (|59|l [or once, if the 
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quartets ((a, /?), (7, 6)) and ((7, S), (a, (3)) are treated as different], and (b) that the central coordinate of the quartet, 
(X°^^,r^ s ), must lie inside the region of the primitive simulation cell. These requirements can be perceived as 
is we have a periodic boundary condition problem with the pairs playing the role of particles and with (X a P ,Y a P) 
serving as the coordinates of the pairs. What can also be learned from expression (|64|l is the fact that k is associated 
with pair-pair correlations. Therefore, its accurate measurement is difficult in systems whose linear L p < 2£, where 
£ is the relevant correlation length. We proceed our discussion assuming that our system is sufficient large and obeys 
the above criterion. 

In order to calculate the third term in Eq. <|59|l we divide our system into an even number of stripes N s = 2M 
(M-integer) parallel to the x axis, as shown in Fig. El The fourth term in Ea. (|59|l is calculated in the same manner 
by dividing the system into the same number of stripes parallel to the y axis. In addition to the primitive cell we 
also need to consider the nearest periodic extensions of linear size L p /2. These periodic extensions, which are also 
shown in Fig. [HJ consist of periodic images of the stripes. We label the stripes included in the primitive cell with the 
numbers M + 1, . . . , 3M, the stripes on the right periodic extension with 1, . . . , M (they are the periodic images of 
stripes 2M + 1, . . . , 3M, and the stripes on the left periodic extension (the images of stripes M + 1, . . . , 2M) with 
3M + 1, . . . ,4M. For each pair we calculate the quantity p afi = <p' (R a P) R^R^ /R af} . The location of the pair, 
which is identified with the mid-coordinate X a/3 = (A Q + A' 3 ) /2, defines the stripe with which the pair should be 
associated. In Fig.|SJeach pair is depicted as a particle. The pair labeled a, for instance, is located in the fifth stripe, 
whereas its periodic image a' is located in stripe number 13. For each stripe i in the primitive cell we calculate the 
sum 

Si = Yl p afS X af3 . (Al) 

pairs in stripe # i 

The sum corresponding to stripe j, the image of stripe is given by 

£,= J2 P aP {X^±L p ), (A2) 

pairs in stripe # i 

where the sign (±) in the above expression depends on whether the image is situated to the right or the left of the 
primitive cell. The product S p S g gives the contribution to the third term in Ea. l|59() of the quartets whose constituent 
pairs are included, respectively, in stripes p and q. These contributions should be in accord with requirements 
(a) and (b), mentioned in the previous paragraph, about the quartets and their locations. In some cases these 
requirements are fulfilled by the image of the stripe rather than the stripe itself. A few illustrative examples are 
given in Fig. |SJ The contribution of the quartets (a, b) and (6, c), for instance, is obtained from the products £5^8 
and SgEn, respectively. The quartet (a, c), on the other hand, should not be introduced into expression 159(1 for 
k via the product S5S11. The distance from a to the image c' is smaller than to c and so the quartet should be 
identified as either (a, c') or as (a', c). The latter is the correct choice because the center of the quartet (a', c) satisfies 

—Lp/2 < X a ' c = (^X a + X c ^j /2 < +L p /2, while the center of the quartet (a, c') falls outside the primitive cell. The 

contribution to the expression for k of this pair is, thus, obtained from the product E11E13. 

The nice feature of the above examples is that the arguments we used to reach our decisions about the correct 
way to handle the quartets have not been based on the precise coordinates of the pairs, but rather on the identity of 
the stripes and their locations with respect to the center of the simulation cell. This means that the products T, p ~E q 
reproduce the contribution of all the quartets corresponding to the relevant stripes. Individual decisions are necessary 
only for a small number of quartets, associated with the following cases: 

• The first case is related with quartets in which the number of stripes separating the pairs is equal to M ', as in 
the case of the pairs b and d in Fig. |SJ which are located, respectively, inside the eighth and the twelfth stripes 
(M = 4 in the above example). The separation between the pairs b and d along the x axis is very close to 
L p /2, and it is impossible to know (without checking the coordinates of the pairs) whether the pair d should 
be replaced by its periodic image d' located in the fourth strips. In a homogeneous system about half of such 
pairs should be exchanged with their images, and so the best estimate for the contribution to expression i|59|) 
for k, arising from quartets including one pair inside the eighth stripe and the other inside the twelfth stripes is: 
0.5E 8 (E 4 + S 12 ). 

• Another case occurs when the stripes containing the two pairs are symmetric with respect to the center of the 
primitive cell and, in addition, the distance between them is larger than M . A typical example is the quartet 
(a, d) in Fig. |SJ in which a is inside the fifth stripe and d is in the twelfth strips. In this case it is obvious 
that (a, d) has to be replaced by either (a,d') or by (a',d), but the two are equally probable. Therefore, the 
contribution of such quartets is is best estimated by: 0.5(S4Es + E12S13) 
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The above rules for correct summation over the different quartets can be summarized by the following compact formula 
for the third term in expression l|59|l : 

I AM AM \ 

EEwO' ( A3 ) 
\p=i 9 =i / 

where the function / is given by 

{1 for \p-q\ < M - 1 and 2M + 1 < p + q < 6M + 1 
0.5 for \p - q\ = M and 2M + 1< p + q < 6M + 1 
0.5 for \p-q\ <M-1 &ndp + q = 2N+l . (A4) 

0.5 for \p - q\ < M - 1 and p + q = 6N + 1 
otherwise 

The value of k obtained using the above expressions [Eas. l|A3l) and l|A4|l ] are not accurate since the contribution of 
some of the quartets is introduced in an approximated way. However, the fraction of such quartets and the resultant 
numerical error can be diminished by taking the limit N s — > oo. In our simulations we have used a set of five 
approximations with N s = 4, 6, 8, 12, 24. 

Another "trick" to speed up the calculation of k: The third and fourth terms in expression (|59|l for k depend on 
the coordinates of the particles. Therefore, several values for these quantities can be obtained from a single MC 
configuration by generating replicas of the original simulation cell. These replicas can be generated by shifting the 
position of the origin of axes, and using the "minimal image convention" to define a replicated primitive cell which 
is centered around the new origin. The computational effort required for the calculation of expression (|59|l in the 
replicas is substantially smaller than that required for the generation of new MC configuration. For one special set of 
replicas the calculation can be done with (almost) no additional effort at all: This set include the replicas generated 
when the origin is shifted by constant intervals Sx = L p /N s in the x direction (Sy = L p /N s in the y direction). Such 
shifts are computationally favorable because they lead to cyclic permutations of the stripes, but do not mix the pairs 
included in each one of them. 
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FIG. 1: A schematic picture of a bilayer membrane (gray) in the reference state (a), and in two deformed states (b) and (c). 
The solid line represents the characteristic surface of the membrane, to which the Helfrich free energy is applied. The areas of 
the characteristic surfaces and the volumes of the membranes (represented by the gray-shaded area in the figure) in (b) and 
(c) are identical. The membrane depicted in (b) is shown in (d) together with the containing cell and the embedding solvent. 
The end points marked by the open circles belong to the perimeter V of the characteristic surface. Another deformation of the 
container, which do not change the total area of the characteristic surface, is shown in (e). 
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FIG. 2: A cylindrical bilayer membrane (gray) with radius of curvature Rq and apex angle 9. The solid line in the middle of 
the membrane represents the characteristic surface. The cylindrical shape of the membrane is obtained via a deformation of 
the containing cell, depicted by the bold dashed line in the figure. The membrane may be thought of as part of a cylindrical 
vesicle (depicted by the thin dashed line) of a similar radius of curvature. 
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FIG. 3: A schematic picture of a lipid molecule in our model system - a trimer consisting of three spherical atoms of diameter 
a. The atom labeled 1 (solid circle) represents the hydrophilic head of the lipid, while the atoms labeled 2 and 3 (open circles) 
represent the hydrophobic tail. 




FIG. 4: Equilibrium configuration of a fluid membrane consisting of 1000 molecules (500 molecules in each monolayer). 




FIG. 5: The inverse of the spectral intensity for undulatory modes l/Z 2 (|/i 9 | 2 ) as a function of the square wave number q 2 . The 
circles mark numerical results, while the solid line depicts Eg, l75H with the values of a and re given by Eg. fro! 
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FIG. 7: The "apparent" bending modulus k* as a function of the inverse of number of stripes dividing the simulation cell, 
1/Ns- The curve depicts the weighted least square fit of a first order (linear) polynomial in 1/N S to the data. 
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FIG. 8: A schematic picture of a system of linear size L v consisting of four pairs (a,b,c,d) and their periodic images (a'.b',c',d'). 
The bold frame marks the boundaries of the primitive simulation cell which is divided into N 3 — 8 stripes labeled from 5 to 
12. The images of the stripes which belong to the nearest periodic extensions of the primitive cell are labeled 1-4 and 13-16 



